Effects of Buckling on Stress and Strain in Thin Randomly Disordered 

Tension-Loaded Sheets 



in 
o 
o 

(N 



(N 



<4-l 

o 

CO 



Bj0rn Skjetne, 1 ' 2 Torbj0rn Helle, 1 and Alex Hansen 2 

'Department of Chemical Engineering, Norwegian University of Science and Technology, 

N-7491 Trondheim, Norway 
^Department of Physics, Norwegian University of Science and Technology, 
N-7491 Trondheim, Norway 
(Dated: February 2, 2008) 

We study how crack buckling affects stress and strain in a thin sheet with random disorder. The 
sheet is modeled as an elastic lattice of beams where each of the beams have individual thresholds 
for breaking. A statistical distribution with an exponential tail towards either weak or strong 
beams is used to generate the thresholds and the magnitude of the disorder can be varied arbitrarily 
between zero and infinity. Applying a uniaxial force couple along the top and bottom rows of the 
lattice, fracture proceeds according to where the ratio of the stress field to the local strength is most 
intense. Since breakdown is initiated from an intact sheet where the first crack appears at random, 
the onset and mode of buckling varies according to where and how the cracks grow. For a wide 
range of disorders the stress-strain relationships for buckling sheets are compared with those for 
non-buckling sheets. The ratio of the buckling to the non-buckling value of the maximum external 
force the system can tolerate before breaking is found to decrease with increasing disorder, as is the 
ratio for the corresponding displacement. 

PACS numbers: 81.40.Jj, 62.20.-x, 05.40.-a 
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I. INTRODUCTION 

In recent years methods have been developed within 
the statistical physics community to describe breakdown 
phenomena in complex media pj . These are the so-called 
lattice models, where the material is reduced to a set of 
points on a grid whereupon disorder is imposed on each of 
the elements on the grid. The desire to understand struc- 
turally non-uniform systems stems from the fact that 
many materials, natural or man-made, show a signifi- 
cant degree of disorder on the microscopic or mesoscopic 
level. In order to realistically describe how such materials 
fracture one has to include the interplay between, on the 
one hand, local variations in material properties and, on 
the other hand, a constantly evolving non-uniform stress 
field. The above mentioned lattice models are especially 
well suited for this purpose. 

Most of the work done with lattice models on frac- 
ture and other breakdown phenomena, however, has fo- 
cused on the fundamental underlying principles rather 
than traditional problems in fracture mechanics. The 
various quantities studied have been expressed through 
scaling laws and critical exponents, often with the aim to 
shed light on universal aspects of phenomena which are 
seemingly unrelated. The most common examples be- 
sides fracture are transport properties and growth pro- 
cesses 0, 0- Obviously there is much to benefit from 
the application of lattice modeling to more specific prob- 
lems in fracture mechanics, especially where disordered 
materials are concerned. 

By far the most popular tool in fundamental studies 
of breakdown processes has been the so-called random 
fuse model a scalar analogue of fracture which re- 
ally models electrical breakdown. Another model, which 



takes account of the vectorial nature of elasticity, is the 
beam lattice Recently, we introduced a three di- 

mensional version of the beam lattice which is suitable 
to describe buckling in thin planar structures 0. Such 
buckling behaviour is perhaps most frequently associated 
with thin plates or beams under compressive loading. In 
this paper we concern ourselves with the special case of a 
thin planar structure under tensile, mode-I type, loading. 
The interaction of buckling with fracture in such cases is 
a well known phenomenon, although as a problem it re- 
mains much less studied. 

Most of the data reported, both theoretical and exper- 
imental, have centered on a few, rather limited, special 
cases, such as that of a thin plate with a center-crack, 
aligned in a perpendicular fashion to the externally ap- 
plied force. When such a plate is subjected to uniax- 
ial tensile loading, transverse compressive stresses build 
up in the vicinity of the crack, causing the unsupported 
edges to deflect out of the initial rest plane. This re- 
distributes the stresses around the crack and leads to a 
stronger singularity at the tip, thus reducing the exter- 
nal force necessary to propagate crack growth. There are 
many applications for which the special case of a homo- 
geneous thin plate with a center-crack is representative. 
Crack buckling, however, is observed under a variety of 
conditions, and often involves anisotropic or disordered 
materials with more than one crack. 

Composites, for instance, are on the increase as the 
preferred material for use in the thin walled plate- or 
shell-structures so essential to the construction of vehi- 
cles for transportation purposes, e.g., hulls and fuselages. 
Critical loads for orthotropic plates have been obtained 
in finite element (FEM) calculations, but only within the 
usual single-crack or hole scenario The importance 
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of buckling and the way it interacts with fracture in the 
presence of multiple cracks has been recognized for some 
time, however. In the aerospace industry, for instance, 
one seeks to make allowance in the design approach for 
the presence of multi-site damage, i.e., assess to what ex- 
tent a series of aligned cracks have on the strength prop- 
erties of a structure 0, 0]. Moreover, buckling plays 
an important role in the breaking of thin sheets where 
the disordered nature of the micro-structure cannot be 
ignored. The crack geometry which obtains in such cases 
may be highly complex. A specific example of this is 
paper. Paper results from a rapid filtration process in- 
volving water and wood fibres. The sheet formed is a 
layered fibre-structure, but nonetheless strongly coupled 
in the vertical direction. Paper is thus a highly stochas- 
tic material, where the essentially random structure is 
modified by flocculation, i.e., an undesired clustering of 
fibres in the early stages of the filtration process 
With conditions of tensile loading frequently arising in 
production facilities as well as printing presses, buckling 
deformations due to tension in paper is a well known 
phenomenon [T^ |. Its interaction with fracture has not 
received sufficient attention, however. 

In the following we briefly summarize some of the 
research that has been done on the buckling of thin 
sheet materials under tensile loading. Among the ear- 
liest investigations was that made by Cherepanov 01 
on membranes containing holes. This, and much of the 
literature which followed mainly concerned itself with 
the calculation of critical loads for the onset of buck- 
ling 0,0,0,0], rather than the effect buckling has on 
the fracture properties once it has set in, i.e., the so-called 
post-buckling behaviour. That buckling should adversely 
affect residual strength has been recognized for some 
time, however, with early experimental observations re- 
ported by Forman [l8|. Dixon and Strannigan |l9j |. and 
Zielsdorff and Carlson 20], hence the interest in deter- 
mining the loads and conditions under which plates with 
specified parameters buckle. As already noted, the source 
to this reduction in strength has been traced to a redis- 
tribution of stresses which leads to a stronger singularity 
at the tip of the crack pHl Ei |23|] . 

Most of the results relevant to the critical buckling 
load have been obtained for thin plates with a center 
crack. Such results are usually expressed in the form 
of an empirical relation which involves plate thickness, 
crack length, Young's modulus and a proportionality fac- 
tor. In their recent experimental work, Guz and Dyshel 
have also considered several cases which can be seen as 
variations on the theme of a central crack; e.g., the ef- 
fect that crack curvature or an inclination angle has on 
either the critical buckling load or the residual strength 
of a plate with a centrally located crack 24]; or the ef- 
fect a straight central crack has on the critical buckling 
load of a two-layered plate j2f| . Centrally cracked plates 
are not the only systems studied, however, plates with 
edge cracks have also been considered. Here the buck- 
ling mechanism has been found to be different from that 



which causes a central crack to bulge |26j . Critical buck- 
ling loads relevant to both perpendicular [27], and 
inclined j2^ edge cracks have been obtained, as well as 
results for the effect buckling has on the residual strength 
of edge cracked panels [23 ■ 

With regard to modeling and theoretical research, an 
early study by Pellet et al. employed a Rayleigh-Ritz 
variational procedure to obtain critical buckling loads in 
the presence of a circular hole [3(J • Recent FEM calcu- 
lations realistically reproduce the observed buckling be- 
haviour around centrally located cracks, and results have 
been obtained for critical loads which agree well with 
experimental findings [H 0SJ |H ^ . The FEM ap- 
proach has also been used to study the various modes 
of buckling and the extent of the buckling zone, e.g., for 
plates with either a perpendicular |2l| or an inclined |2^] 
central crack. Gilabert et al. also obtained results rele- 
vant to the zone of deformation [3l| , and critical loads for 
various crack geometries, e.g., circular holes or rectangu- 
lar cut-outs with sharp or rounded corners, have been 
obtained in other FEM calculations [32| . 

Features of the post-buckling behaviour, other than 
the shape and extent of the buckling zone, was ob- 
tained by Petyt, i.e., for the vibration characteristics of 
a centrally cracked plate subject to acoustic loads |21j . 
Petyt also addressed the non-linear nature of FEM cal- 
culations for the post-buckling behaviour, and Riks et 
al. |22j used such an analysis to show that the energy 
release-rate at the tip of the crack undergoes a sudden 
increase at the onset of buckling. The stress intensity 
along the post-buckling path is then larger than that ob- 
tained along the pre-buckling path for the same load, a 
result which, in agreement with experimental observa- 
tions, indicates that the residual strength of the plate is 
reduced by buckling. The effect of crack inclination on 
the energy release-rate in the post-buckling state has also 
been studied j^. FEM calculations for the load versus 
crack-opening length in buckling and non-buckling frac- 
ture modes have been carried out in a study by Seshadri 
and Newman showing a significant reduction in the 
residual strength. Their work also considered the effect 
of plasticity by assuming a hypothetical material with a 
very high crack-tip opening angle, with the reduction in 
strength due to buckling now being generally less pro- 
nounced than in the brittle case. 

As the above summary shows, practically all previous 
work considers the effect buckling has on the strength 
properties of an already cracked plate, or a plate with a 
geometrical discontinuity such as a circular hole or a rect- 
angular cut-out. In other words, if the physical parame- 
ters of the plate are such that buckling can be expected 
before the crack begins to grow, the residual strength of 
the plate will be significantly lower than what would oth- 
erwise be expected, based on an analysis which does not 
take account of buckling. The present study of fracture 
and buckling is fundamentally different in the sense that 
we regard a sheet which, in its initial state, has no cracks 
or other discontinuities. Instead, cracks form by a com- 
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plex process which depends on the evolving distribution 
of stresses and its interaction with a disordered meso- 
structure. The onset of buckling in this scenario, and 
the effect buckling has on the fracture properties, will 
vary according to the type of disorder used, i.e., weak or 
strong. Whereas for strong disorders there will be signif- 
icant sample-to-sample variations, such variations tend 
to be less pronounced for weak disorders. However, even 
for weak disorders the final crack which breaks the sys- 
tem will only rarely appear at the exact center of the 
sheet, and even then the situation might be complicated 
by additional cracks in the vicinity - cracks which inter- 
act with the main crack so as to alter the distribution 
of stresses and hence also the exact shape or mode of 
buckling. Therefore, due to the statistical nature of the 
results obtained, features such as the extent of the buck- 
ling zone, or the shape of the deflected crack edge, will 
not at present be dealt with in any detail. For the same 
reasons critical loads are not calculated, since the mag- 
nitude of this quantity depends on very specific sheet 
parameters, i.e., for a given sheet thickness the critical 
load has been shown to depend on the ratio of the crack 
length to the sheet width. 



II. THE BEAM LATTICE 

The beam lattice used in our calculations is a reg- 
ular square lattice, where each beam has unit length. 
System size L therefore corresponds to the number of 
beams along the top or bottom rows. The nodes are 
equidistantly spaced along J = L + 2 horizontal rows 
and I = L + 1 vertical columns, each having four near- 
est neighbours to which it is fastened by elastic beams. 
When nodes are displaced the angle at the joint where 
two beams come together remains perpendicular, thus in- 
ducing shearing forces and bending moments in addition 
to axial tension or compression. 

In the plane beam lattice there are three degrees of 
freedom for the displacement of nodes, i.e., translation 
along either the X-axis or the Y-axis, and rotations 
about the Z-axis. The displacement field is obtained by 
requiring the sum of forces and moments on each node 
to be zero. Specifically, we solve 
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by numerical relaxation, i.e., the conjugate gradient 
method [3^, to obtain the set of displacements which 
minimizes the elastic energy of the lattice. 



In Eqs and 10 , A and T denote axial and transverse 
force, respectively, while in Eq. 10} M denotes the bend- 
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ing moment. Hence, X A\ is the force exerted on node i 
from j '■ — 3 along the X-axis by axial tension or compres- 
sion. Neighbouring nodes are numbered anti-clockwisc, 
starting with j = I on the left. 

Defining Sr — rj — n, where r £ {x,y,w}, the contri- 
butions from j — 1 are 



x A^ = Ux, 



[Sy- -(Wi + Wj)]. 
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are the prefactors characteristic of the material and its 
dimensions, i.e., E is Young's modulus, p and I the area 
of the beam section and its moment of inertia, respec- 
tively, and G the shear modulus [fj. 

The fracture process consists of removing one beam at 
a time, whereby a new set of displacements are obtained 
at each step by solving Eq. Q. The criterion by which a 
beam is removed from the lattice depends on the ratio of 
the local stress to the breaking threshold. Using and 
tjvf for the maximum thresholds in axial force and bend- 
ing moment, respectively, a good breaking criterion [(| 
inspired from Tresca's theory is 



A_Y + \M\ 
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where \M\ = maxflM,- |, |Mj' ; |) is the largest of the mo- 
ments at the two beam ends i and j. 

The time taken for mechanical equilibrium to be 
reached is assumed to be much shorter than the time 
taken to remove a beam, i.e., the fracture process is as- 
sumed to be quasi-static. It is driven by imposing a fixed 
unit displacement on the top row of the lattice. Since 
internal displacements, forces and moments are propor- 
tional to this, the actual external elongation of the lat- 
tice is obtained by determining the minimum value of the 
proportionality constant Al in 
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from which the external force is obtained as 
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FIG. 1: A disordered lattice of size L — 20, shown at four dif- 
ferent stages in the breakdown process. The lattice is strained 
to failure by applying a force couple at the top and bottom. 
With the appearance of large cracks, the structure is seen to 
deflect out of the initial rest plane. The number of broken 
beams in this simulation are, from (a) to (d), N = 65, 72, 78, 
and 86, respectively. 

keeping track of whether beams are intact (1) or have 
been broken (0), respectively. 

In Eq. Hill) , contributions other than y A^ cancel when 
the sum is over the entire lattice. This is due to the 
square lattice topology and the nature of the external 
boundary conditions applied, i.e., mode-I type loading in 
the ^-direction. Regarding the internal forces, the same 
term continues to be the sole non-zero contribution when 
the lattice is intact. Consequently the first beam to break 
is that for which the ratio A/tA is the largest. After this 
has been removed, however, bending moments M and 



transverse forces T (including shear) are induced in the 
immediate neighbourhood of the beam. This is due to 
the screening effect of the "hole" , or crack, created by its 
removal from the lattice. 

A thin sheet will usually display deviations in symme- 
try with respect to the thickness, e.g., there may be vari- 
ations in the thickness itself or there may be a gradient 
in the structural properties of the material. An example 
of the latter is paper, where, due to the process by which 
it is manufactured, the fibre structure on one side always 
has a stronger orientational bias. In other materials the 
density varies in the thickness direction. When a uniax- 
ial force couple is applied on opposite edges of the sheet, 
such variations create bending moments about axes in 
parallel within the ATF-plane, see Fig. ^ which shows 
the coordinate system and the direction of the external 
load. In fact, when the internal stresses (which arise as 
a consequence of the external load condition) combine 
with certain crack configurations, minute deviations of 
the symmetry plane itself from a perfect two dimensional 
embedding will be sufficient to cause buckling. Numerous 
studies have been reported in the literature concerning 
the external load necessary to cause buckling, e.g., for a 
sheet with a central crack the magnitude of the critical 
load has been found to decrease with decreasing sheet 
thickness and increasing crack extent. 

The additional terms which cause buckling are much 
smaller in magnitude than those governing the forces 
within the plane lattice. There is, however, a non- 
separable relationship between in-plane and out-of-plane 
displacements which causes the in-plane coordinates of 
the non-buckling lattice to change significantly when 
buckling is allowed. For this reason the X i: and Wi 
components of the buckling lattice Q 





' Ui ' 




'Ui - 




Vi 




Vi 




Wi 
Xi 


= X 


w t 


3 


Vi 










Y 








.Zi _ 



contain additional non- linear terms, i.e., terms not in- 
cluded in Eq. (J). 

Specifically, the axial force component 

Xp =-Sx (14) 
a 

in Eq. (JIJ is replaced by 
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where 



(16) 



is the force along the axis of the beam, including an angular correction which takes into account the additional 
elongation due to bending. Likewise, the transverse force 
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for the displacements normal to the XY-plane, and 
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for the rotations about the Y-axis. Finally, 
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is the torque of the beam when rotations are about the 
X-axis. Assuming w > t, 



(24) 



is the torsional moment of inertia in Eq. (|23[) . with w 
denoting the width of the beam cross section and t its 
thickness. Assuming a rectangular cross-section, the mo- 
ments of inertia for bending are 
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within the AK-plane, and 
1 



(25) 



(26) 



within the YZ- and XZ-planes, respectively. 

The expressions for the forces acting on the beams in 
Eqs. (|15fl to i|23|) have been derived by considering an 



elastic beam with no end restraints |34| . where the ratio 
of the beam width to the thickness presently has been set 
to 10:1. With regard to bending flexibility, the lattice is 
now more pliable in the out-of-plane direction, as would 
be expected for a thin sheet material. 

In lattice modeling the rule by which a beam is bro- 
ken can be specified according to the properties of the 
material one wishes to study. Presently the fracture cri- 
terion is taken to depend on a combination of axial stress, 
bending and torsion. Hence, we assume 
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is the effective stress, and 
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is the torque. Moreover, with /_tc denoting the combined 
bending moment and 
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being the aspect ratio of the cross section of the beam, 
the expression 
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(31) 



is an enhancement factor in Eq. I|28|) . In Eq. 129|) the 
Kronecker delta has been used to distinguish between 
the four neighbouring beams. 

Angular displacements about the X- or F-axis in 
Eq. (|27|l activate the stress enhancement mechanism. 
This increases the stress in a beam when it is under axial 
tensile loading. Specifically, the larger the the load is, 
the more sensitive the beam will be to the presence of 
a certain amount of axial torque. Compressive loads are 
assumed to be less important, with torque now instead 
removing some of the axial compression. 

The exact mechanism by which buckling alters the rup- 
ture mode of a thin sheet will probably vary according to 
material properties, structural composition, and so on. 
In many cases it is a well known fact that the work re- 
quired to drive a crack across a given area is much smaller 
in mode-Ill tearing than in pure mode-I tension, as is 
easily verified with a piece of paper. Transverse forces, 
however, are not presently assumed to contribute. This is 
because the disorder of the sheet is modeled on a meso- 
scopic scale. In a material such as paper, tearing is a 
shear displacement which affects material properties on 
much smaller scales, e.g., on the level of the individual 
fibres. In the beam lattice, on the other hand, each in- 
dividual beam is representative of the sheet on the level 
of fibre flocculations. The effect of mode-Ill crack prop- 
agation is instead included by the above combination of 
torsion and axial stress. 



Since the leading terms of Xj, Yi and Wi are all linear, 
the actual solution for this in-plane projection always lies 
close to its linear solution. It is found by re-initializing 
the search, each time using conjugate gradients starting 
from the previous linear approximation. This is repeated 
until the minimum stops changing, typically 6-7 searches 
are required, with convergence being rapid. 

After this intermediate solution has been obtained, it 
is frozen, whereupon the out-of-plane coordinates are re- 
laxed, one at a time. In this case, however, the lead- 
ing terms are non-linear. Consequently a single search is 
made toward the minimum to obtain 
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i.e., a partial only, or (at best) very approximate, solu- 
tion. Moreover, in order to ensure that this incomplete 
move carries towards (and not away from) the minimum, 
the step size of the conjugate gradient iterations in this 
phase is reduced to a much smaller value. 

Likewise, the out-of-plane angular displacements are 
updated, one at a time, using the same down-scaled step 
size, to obtain 



00 



00 



(34) 



while keeping all other coordinates fixed. 

After re-setting the itcrational step size, the whole pro- 
cedure outlined above is repeated. The updated coordi- 
nates obtained for the out-of-plane displacements, ap- 
proximate as they are, do nonetheless cause the in-plane 
displacements to change. As the final buckled config- 
uration of the lattice is approached the quality of the 
intermediate partial solutions, represented by Eqs. (|33|) 
and <|51|l . gradually improves. 

Hence, after a number of repetitions we obtain 



III. NUMERICAL SCHEME OF 
CALCULATIONS 

Mathematically, conjugate gradients is an iterative 
method to obtain the minimum of a quadratic expres- 
sion, in our case the elastic energy. For the energy to be 
quadratic, however, the forces involved must be linear. In 
obtaining a numerical solution, therefore, the presence of 
non-linear terms is a complicating factor. Nonetheless, 
provided the proper numerical safeguards are employed, 
the correct minimum can be found effectively by use of 
conjugate gradients. Specifically, 
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is the solution obtained by relaxing the in-plane coor- 
dinates while keeping the out-of-plane coordinates fixed. 
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for the sum of forces and moments on all nodes. The pre- 
vious set of displacements is now identical to the current 
set of displacements, the calculation having converged 
upon the final solution. 



IV. DISORDER 

Each time a beam is broken, a new set of displace- 
ments is calculated according to the scheme outlined in 
section lnTl A fundamental factor deciding how the lattice 
breaks is the choice made for the type and magnitude of 
disorder in the distribution of breaking thresholds. One 
of the reasons why lattice models are practical is the ease 
with which such disorder may be included. 
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FIG. 2: The lower remaining part of a beam lattice of size L = 50 after it has been broken completely, shown for five different 
disorders, i.e., from left to right: D = 0.25, D — 0.5, D = 1, D = 2 and D = 4, respectively. From top to bottom, six different 
samples have been included in each case, the only difference between the samples being the random casts generated for the 
breaking thresholds. 



Presently we generate a random number r on the unit 
interval [0, 1] and let this represent the cumulative thresh- 
old distribution. In Eq. I|27|l the breaking thresholds are 
now assigned as t = r D , with 

Pit) = (36) 

being the probability density. The same distribution is 
assumed for the threshold in axial force, t = tp c , and 
bending moment, t = i M c> with the random casts, how- 
ever, being different in the two cases. 

There are now two types of distribution, i.e., D > 0, 
in which case 

P(t) - fh (37) 



is a cumulative distribution with bounds < t < 1, and 
D < 0, in which case 

P{t) = l-ti (38) 

is a cumulative distribution with bounds 1 < t < oo. 
In this prescription D = corresponds to no disorder. 
As |D| increases the coefficient of variation with respect 
to any two random numbers r and r' on the interval [0, 1] 
also increases, with the coefficients for D > and D < 
being reciprocal but otherwise the same. Hence, large 
values of \D\ correspond to strong disorders and small 
values to weak disorder. A few examples for D > have 
been included in Fig. [3 where the bottom part of the 
broken lattice is shown for five different disorders. Also 
shown is the sample-to-sample variation for each of the 
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A 



FIG. 3: Force / versus displacement A for an L = 32 lat- 
tice with disorder D — 1. The two upper curves denote 
non-buckling fracture, with the continuous line representing 
Eq. Q and open circles representing Eq. 11311 with buckling 
suppressed. The lower curve (filled circles) denotes buckling 
fracture, i.e., the results of Eq. 1131 with the out-of-plane 
degrees of freedom now included. 

disorders. For a given magnitude of D, the position of 
the final crack which breaks the system is seen to vary 
from one sample to the next, as is its morphology. As 
the disorder increases, so does the roughness of the crack 
interface. The number of beams removed also increases 
with the disorder, i.e., the D — 4 samples are seen to 
be somewhat more diluted than the D = 2 samples. In 
Eq. i|3(j|) . D > and D < represent widely different 
types of distribution. While the former is a power law 
with a maximum threshold of one, and a tail which ex- 
tends toward zero, the latter is a power law bounded be- 
low by a minimum threshold of one, but now with a tail 
which extends toward infinity. Both D > and D < 
are included in the present calculations. 

In the past many different distributions have been used 
to generate random breaking thresholds. However, as 
shown by Hansen et al. |35|. as the system size diverges 
only the power law tails of the distribution, if they ex- 
ist, towards zero or infinity should matter. Hence the 
use of D as a parameter is very convenient, enabling the 
asymptotic behaviour of the fracture process to be fully 
explored as a function of the disorder. 



V. STRESS AND STRAIN 

In the absence of structural disorder the crack now 
grows laterally from the site of the first beam removed, 
taking the shortest possible path across the lattice. Since 
in our model the beams are linearly elastic up to the 
breaking threshold, the first break triggers catastrophic 
rupture. Stress and strain evolves differently in the pres- 
ence of disorder. Now there are two competing mecha- 
nisms for crack growth. On the one hand, the presence 
of a crack causes stress to be intensified in its immediate 
vicinity, thereby lending bias towards the growth of al- 



ready existing cracks. On the other hand, variations in 
material strength dictate that new cracks should instead 
appear in regions which are structurally weak. Which of 
the two mechanisms is the most important depends on 
the disorder regime. While in the case of strong disorder 
fracture is initially disorder dominated, it tends to be lo- 
calized from the very beginning in the case of weak disor- 
der. For strong disorders small cracks appear at random 
in the early stages of the process. Here the dominat- 
ing feature is a wide distribution of breaking thresholds. 
Since the weakest beams tend to be removed first, the dis- 
tribution gradually narrows as more beams are removed. 
Simultaneously, with a growing number of cracks appear- 
ing on the lattice, a highly non-uniform stress field devel- 
ops. In other words, the distribution of stresses widens. 
At the point where the fracture process goes from being 
disorder dominated to stress dominated, crack growth 
becomes localized |35j . Smaller cracks now merge into a 
single dominating crack and the evolution of stress with 
strain goes from being stable to unstable. 

For a system of size L = 32 and disorder D = 1, a com- 
parison between the buckling and non-buckling stress- 
strain characteristics is shown in Fig. |3J The average 
stress and the average strain has been computed for ev- 
ery beam broken, and the number of samples involved is 
10000 in the non-buckling case and 975 in the buckling 
case. Also included is the result of Eq. (|13(l with the 
out-of-plane degrees of freedom suppressed. This result 
is based on 525 numerical realizations. The agreement 
between the non-buckling results of Eqs. and (|13|) 
is seen to be excellent, especially in the controlled and 
early catastrophic regimes. Towards the end of the catas- 
trophic regime the loads obtained with Eq. I|13|l are very 
slightly lower than those that are obtained with Eq. JQ, 
a result which can be ascribed to the presence of non- 
linear terms in the former. With buckling, a significant 
reduction is obtained in both maximum strength and 
displacement. There is also a notable difference in the 
shape of the curve within the catastrophic regime. Here 



TABLE I: Ratio of buckling to non-buckling maxima, ob- 
tained for the external displacement A and force /, for disor- 
der D = 1. The total number of samples calculated is S, and 
L is the system size. 



L 


Az/Ao a 


hlh 


S z 


So 


14 


0.83 


0.93 


1500 


5000 


17 


0.92 


0.94 


500 


2500 


20 


0.79 


0.92 


1000 


1000 


23 


0.83 


0.92 


203 


800 


27 


0.83 


0.91 


328 


600 


32 


0.75 


0.90 


975 


10000 


40 


0.74 


0.89 


210 


1400 


50 


0.77 


0.91 


70 


1750 


63 


0.77 


0.89 


110 


700 


80 


0.80 


0.92 


55 


550 



"Quantities labeled Z refer to the buckling case. 



9 



the response is less stable with respect to displacement 
control. That is, the force in the catastrophic regime falls 
off more rapidly as the displacement increases. 

Stress and strain for a range of system sizes is shown 
in Fig. 0] for the same disorder, i.e., D = 1. In calcula- 
tions for the non-buckling beam lattice, involving a much 
larger range of sizes j^, the scaling with L of the top 
of the stress-strain curve is found to be characterized by 
an exponent close to unity. The stress-strain curves can 
then be made to collapse onto each other by scaling the 
axes according to 



f/Ui = ^(A/L 7 ), 



(39) 



where 7 « 1 and cj> is a scaling function. Since there 
is no reason why the buckling system should behave ac- 
cording to different laws in this respect, the reduction 
in stress and strain should itself be proportional to sys- 
tem size. As was noted in section |n] fracture is initiated 
by imposing on the top row of the lattice a displace- 
ment of one beam length. Hence, to avoid scale effects 
on the buckling behaviour, one of the factors L intro- 
duced in the stress enhancement factor, i.e., in Eq. (|31JI . 
is a scale factor. Without this factor, a different value 
of the exponent 7 would be obtained in Eq. (|39[) . With 
the current choice of parameters, maximum stress and 
strain in the buckling and non-buckling cases scale ac- 
cording to the same law, as can be seen from Fig. 0] and 
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FIG. 4: Force / versus displacement A for a range of system 
sizes with disorder D = 1, i.e., for (a) L = 23, (b) L = 27, 
(c) L = 32, (d) L = 40, (e) L = 50 and (f) L = 63. In each 
case the top curve is the non-buckling result of the simple 
beam model, calculated from Eq. and the curve below 
is the buckling result, calculated from Eq. 1131 . The labels 
on the axes are scaled down from those in plot (f), being 
otherwise proportional to system size. 



/ 




FIG. 5: Force / versus displacement A for an L — 32 lattice 
with disorder D = 1, comparing the buckling (filled circles) 
and the non-buckling (open circles) results of Eq. il l-'-il where, 
in the former case, x = in Eq. 1281 . 



also from the comparison of buckling and non-buckling 
stress-strain maxima in Table [I] Here the values obtained 
for the reduction in maximum stress and strain appear 
to be consistent for systems larger than about L = 20. 
Below this, finite size effects become apparent. The most 
reliable estimates are obtained with the largest number 
of calculated samples, hence, for D = 1, buckling reduces 
the maximum strength by about 10% and the maximum 
displacement by about 25%. The shape of the curve in 
the catastrophic regime varies according to system size. 
Beyond the turn-over point between stable and unstable 
crack growth, rupture in the non-buckling system is seen 
to become increasingly less stable as the size of the sys- 
tem increases. This effect is even more pronounced when 
the sheet is allowed to buckle. 

Significant differences are evident in a comparison be- 
tween the force or displacement fields of, say, a uniform 
center-cracked lattice in the case of buckling with the 
corresponding force or displacement fields in the non- 
buckling case. Specifically, the transverse forces near 
to the crack edges, which are compressive in the non- 
buckling lattice, are released when the lattice buckles, 
causing the flanks of the crack to deflect. Since the alter- 
ations in the force or displacement fields extend beyond 
the immediate neighbourhood of the crack tips, one may 
ask whether these effects in themselves are sufficient to 
bring about a reduction in the maximum load carrying 
capacity of the lattice. The mechanism by which the 
stress is intensified at the crack tips, however, takes place 
on a scale smaller than the individual beam, which is 
why the fracture criterion Eq. I|27|) has been augmented 
by the factor x- Hence, the hypothetical case of fracture 
where buckling does not induce intensified stress at the 
crack tips can be investigated simply by setting x = in 
Eq. (|28l) . In Fig. [S] the result is compared with that ob- 
tained in non-buckling fracture. Clearly, the evolution of 
stress with strain is seen to be the same in both the sta- 
ble and catastrophic regimes. Both curves were obtained 
from Eq. (|13fl . based on 1350 samples in the buckling case 
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0.1 




TABLE II: Ratio of buckling to non-buckling maxima, ob- 
tained for the displacement A and force /, for L = 32. The 
number of samples is 5*, and D is the disorder. 



20 40 60 80 100 1.8 3.6 5.4 7.2 9.0 
A 

FIG. 6: Force / versus displacement A obtained for a lattice 
system of size L = 32, for a range of disorders with both 
D < and D > 0. On the left are shown (a) D = -5, 
(b) D = -3, (c) D = -2 and (d) D = -1. On the right are 
shown (e) D = 0.33, (f) D = 0.5, (g) D = 1 and (h) D = 2. 
In (a)-(c) the respective scales on A are x5, x2.5 and xl.5 
that in (d), and in (e)-(g) the respective ratios //A are 0.38, 
0.3, and 0.15. In all cases the stress-strain curve for buckling 
lies below that which does not buckle. 



and 525 samples in the non-buckling case. One can thus 
state with certainty that it is the intensified stress at the 
crack tips, due to a coupling between in-plane and out- 
of-planc deformations, rather than the re-distribution of 
stresses within the buckling zone (but away from the im- 
mediate neighbourhood of the crack tips) which causes 
the reduction in residual strength. 

In Fig is shown the effect disorder has on the inter- 
action of buckling with fracture. Plots on the left-hand 
side display stress-strain curves for D < type disor- 
der. Here the initial response, i.e., the linear relationship 
which extends from the origin to the data point of the 
first broken beam, has been omitted. There is no tail 
towards zero in the distribution of thresholds here, and 
consequently there are no broken beams on this part of 
the curve. In plot (d), where D = — 1, the first beam to 
break triggers a catastrophic rupture mode which back- 
tracks along the initial linear response for the first few 
breaks. It then encounters a vertical section of the curve 
where several values of / correspond to the same A. If 
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10000 
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192 


1000 



"Quantities labeled Z refer to the buckling case. 



displacement control is applied and relaxed sufficiently 
fast for the crack to be halted at this point, we have a 
situation of conditional stability where a slight perturba- 
tion, say a bump or a jar, suffices to further propagate 
the crack (this refers to the average situation, with in- 
dividual samples being subject to fluctuations). Other- 
wise, applying displacement control without this sudden 
relaxation, the crack develops catastrophically until it is 
arrested when encountering strong beams in the tail to- 
wards infinity. From here on, the force continues to fall 
off as the displacement is increased. 

The main effect buckling has for weak D < disorder is 
to make the force fall off more rapidly in the catastrophic 
regime. Additionally, the section of the curve which is 
conditionally stable in the non-buckling case is rendered 
unstable, i.e., the curve turns back on itself. As |D| in- 
creases there is a turn-over in the average stress-strain 
behaviour in the sense that, beyond D — —3, force con- 
trol may be applied without necessarily triggering catas- 
trophic rupture. This, of course, is due to the presence 
of a large number of beams with high breaking thresh- 
olds. When the tail towards infinity becomes sufficiently 
strong, in other words, the number of beams which can 
be found in the vicinity of the lower bound becomes a 
minority. The stress-strain relationship then attains a 
similar form to that of D > 0, except now fracture starts 
at a finite displacement or force. Although in (a) the con- 
trolled regime, which obtains after the first beam breaks, 
contains a smaller number of broken beams than does, 
for instance, the one in (g), the reduction in force and 
displacement due to buckling is comparable in the two 
cases. The reason for this is a more intense stress field 
in the former case, caused by higher thresholds, which in 
turn moves the onset of buckling to an earlier stage of 
the fracture process. 

Displayed on the right-hand side are stress-strain 
curves with D > type disorder. These are mostly sub- 
ject to the same features as the result of Fig. |3| relevant 
to D = 1. An exception, perhaps, is D = 2, for which 
the stability in the catastrophic regime appears to be 
unchanged by buckling. For D = 2 and beyond, how- 
ever, the number of beams relevant to the catastrophic 
regime is small compared to that of the controlled regime. 
This means that the number of samples which contribute 
decreases toward the end of the stress-strain curve (the 
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curves have not been truncated at the average number of 
broken beams), and hence statistical fluctuations become 
large in this region. 

Whereas for weak D > disorders only a small re- 
duction is obtained in the maximum of stress and strain, 
the stability in the catastrophic regime of fracture is sig- 
nificantly affected in this case, as can be seen from (e) 
in Fig. i.e., for D — 0.33. The reason is the onset of 
buckling, which for low disorders occurs near the top of 
the curve. Even when the disorder is sufficiently low for 
the onset of buckling, in average, to occur after the top 
has been reached, a slight decrease in maximum strength 
may be expected. This, of course, is due to the fact that a 
number of samples will buckle prior to this average onset. 

In Table [H] results for the L = 32 system are shown 
for a range of disorders with D > 0. Here the decrease 
in force and displacement is seen to depend on the mag- 
nitude of D, i.e., as \D\ increases buckling has an in- 
creasingly adverse effect on both the maximum load and 
the maximum displacement a disordered system can sus- 
tain. The maximum displacement is more strongly af- 
fected than the maximum load. 



VI. SUMMARY 

The breaking characteristics of thin sheets with struc- 
tural disorder have been obtained in numerical simula- 
tions which include the out-of-plane buckling behaviour. 



The model used is an elastic lattice of beams where each 
beam is representative of the scale of the structural disor- 
der. Depending on the magnitude of disorder, breakdown 
is either localized to the first point of damage or initially 
a random cracking process which at a later stage crosses 
over to localized fracture behaviour. 

The breakdown process is initiated from an initially in- 
tact sheet, where buckling sets in after a certain amount 
of damage has occurred. Specifically, the onset of buck- 
ling varies considerably according to both the size and 
configuration of the emerging cracks. Given a certain 
system size and disorder, several numerical realizations 
of a sheet are generated, corresponding to different sets 
of random breaking thresholds. The statistical proper- 
ties are then obtained from the average behaviour based 
on the disorders and sizes chosen. 

As in the case of uniform pre-cracked sheets, it is found 
that buckling adversely affects the external force and dis- 
placement a randomly disordered sheet can sustain in 
mode-I type tensile loading. The degree to which the 
maximum force and displacement is reduced depends on 
the magnitude of the disorder. For instance, in a material 
such as paper this would mean that buckling should af- 
fect the maximum load carrying capacity more adversely 
in the case of a fibre-web with uneven formation than 
one with a more even formation. When the meso-scale 
disorder is low the reduction in strength is insignificant 
and it is the catastrophic regime which is most affected, 
now being less stable. 
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